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The evolution of infinitesimal, localized perturbations is investigated in a one-dimensional di- 
atomic gas of hard-point particles (HPG) and thereby connected to energy diffusion. As a result, 
a Levy walk description, which was so far invoked to explain anomalous heat conductivity in the 
context of non-interacting particles is here shown to extend to the general case of truly many-body 
systems. Our approach does not only provide a firm evidence that energy diffusion is anomalous 
in the HPG, but proves definitely superior to direct methods for estimating the divergence rate of 
heat conductivity which turns out to be 0.333 ± 0.004, in perfect agreement with the dynamical 
renormalization-group prediction (1/3). 
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After the discovery of anomalous heat conductivity in 
classical one-dimensional lattice systems P|, in the last 
years a renewed attention has been devoted to the old 
problem of identifying the minimal ingredients required 
for the Fourier law to be ensured. As summarized in 
a recent review article 0, many different models have 
been numerically investigated to identify the physical 
conditions under which the thermal conductivity k di- 
verges with the system size L and, having assessed that 
K w L", to determine the possibly different universality 
classes for the divergence rate a. Simultaneously, sev- 
eral attempts have been made to estimate analytically 
the scaling behaviour of n: self-consistent mode-coupling 
theory[3| and the Boltzmann equation suggest that 
a = 2/5, while dynamical renormalization group indi- 
cates a = 1/3 4]. Both predictions are compatible with 
numerical simulations which are, however, often affected 
by relatively strong finite-size corrections. The only sys- 
tem where convincing results have been obtained is the 
Fermi, Pasta Ulam /3-model in the infinite temperature 
limit. Its behaviour is consistent with a — 2/5 5], but 
the simmetry of the potential casts doubts about the gen- 
erality of this model|3]- A further simple system that can 
be effectively simulated on a computer is the diatomic 
hard-point gas: there, interactions are provided by elas- 
tic collisions of point-like particles'^. Unfortunately, the 
most detailed numerical simulations reported in the lit- 
erature show a slow growth of the divergence rate with 
L, so that some authors claim that a = 1/3 j^, while 
others state that the conservative guess a = 0.25 is more 
realistic 0- Settling this issue is not only conceptually 
important, but it is a necessary requisite to later quan- 
tify finite-size corrections, a crucial issue in applications 
to, e.g., carbon nanotubes, where one needs to know the 
prefactor as well. 

Although the problem involves intrinsically many de- 
grees of freedom, some researchers have tried to shed 



some light with reference to the simpler setup of non- 
interacting particles moving along a periodic array of con- 
vex scatterers (billiard gas channels) The absence 
of interactions simplifies the task of understanding heat 
conductivity and allows, in particular, tracing back heat 
conductivity properties to the diffusion of single particles 
at equilibrium. Assuming that the mean square displace- 
ment {x'^{t)) scales as {(3 — 1 corresponds to normal 
diffusion), it can be shown that a = (3 — 1 0, under 
the assumption that each particle exchanges energy only 
at the channel borders, where thermal reservoirs operate. 
The limit of this approach is, on the one hand, that the 
relationship between the diffusion exponent (3 and the mi- 
croscopic dynamics remains to be established [l^ and, on 
the other hand, that particles do not mutually interact. 
Nevertheless, in this Letter we show that upon interpret- 
ing the energy density as a pseudo-particle density, the 
above reasoning can be fruitfully applied to the HPG to 
find a convincing evidence that energy diffuses like in a 
Levy walk process |l^. thus bridging two research lines 
that were so far basically disconnected from one another. 
As a further result, we are also able to establish that in 
the HPG, a = 1/3 with a 1% accuracy. 

The model consists in a chain of N point-like parti- 
cles with alternating masses mi lying on a segment of 
length L, Q < Xi < L (i = \, . . . ,N), where m2i — m, 
with (i = 1, . . . , [A^/2]), and 17121+1 = rm, with (i = 
0, . . . , [(TV — l)/2]), where square brackets indicate the 
integer part. Due to the absence of intrinsic energy- and 
length-scales, we can fix them at wish. The only relevant 
parameter that cannot be scaled out is the mass ratio 
r. Without loss of generality, we set the number density 
V = N/L — 1, the energy per particle e = {mivf)/2 = 1 
(which, in turn, fixes the temperature, ksT = 2) and 
one of the two particle masses equal to 1. The dynamics 
of this model is very simple, since the velocities change 
only in the collisions between adjacent particles while the 



updating rule is determined by the conservation laws of 
kinetic energy and linear momentum. By denoting with 
Vi (vl) the velocity before (after) a collision, the evolution 
equation amounts to 

1 ~ r 

= «j ± "^j) (1) 

where a plus (minus) sign is selected depending whether 
i is even (odd) and j = i ± 1. It is also interesting to 
notice that the positions Xi contribute only indirectly to 
the evolution, by determining the collision times. Such 
properties, plus the conservation of the particle ordering 
along the chain, allow simulating the dynamics with an 
event driven algorithm that exploits the heap struc- 
ture of the future collision times. 

When r = 1, collisions just exchange particle velocities, 
that are, thereby, integrals of motion. Away from r — 
1, the system is no longer integrable, but remains non- 
chaotic. In fact, since the evolution equation is linear, 
it describes the dynamics of an infinitesimal perturbation 
dvi{t) as well. Therefore, the weighted Euclidean norm 

Q = E'^(2)(''^) = E"^*('^«*w)' (2) 

i i 

of a generic perturbation in tangent space is conserved 
for the same reason that energy is conserved in the phase- 
space dynamics, and, a fortiori the maximum Lyapunov 
exponent cannot be larger than zero 

A more detailed characterization of the dynamics can 
be obtained from the spectrum A{u) of convective Lya- 
punov exponents [l^ . which quantifies the space-time 
growth rate of a perturbation i5(2)(*,i = 0) = Si^ ini- 
tially localized in i = {Sij is Kronecker S function) 

A(u) = lim i log 5(2) [i = ut, t) . (3) 

A{u) attains its maximum at u = 0, where it coincides 
with the maximum Lyapunov exponent. In a standard 
chaotic model, upon increasing A(m) decreases until 
it crosses at m = u^, a value that has been shown to co- 
incide with the sound velocity 01 • then interesting 
to understand how the scenario modifies in a non-chaotic 
model such as the HPG. Previous studies revealed a slow 
convergence to the asymptotic regime both when r is 
close to 1 and for r 1 |3| ■ In order to avoid such prob- 
lems we have chosen r = 3, but we can stress that the 
same scenario has been confirmed by a few simulations 
run for r = 2. The data reported in Fig. 1 shows that A 
converges to for u < Uc ~ 1.02 (uc has been determined 
by fitting the long-time behaviour of the secondary maxi- 
mum of A(u)), while it becomes strictly negative at larger 
velocities. This is analogous to harmonic chains, where 
the exponential growth rate of a perturbation is strictly 
zero as long \u\ is smaller than the sound velocity [15|. 
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FIG. 1: Spectrum of comoving Lyapunov exponents at differ- 
ent times for r = 3; the chain length is N=8190. From bottom 
to top the curves have been obtained by comparing pertur- 
bations ad times (100,200), (200,400), (400,800), (800,1200) 
and (800,1600). The perturbation amplitude has been aver- 
aged over 10* different realizations. 

We have verified that Uc coincides, within the numerical 
error, with the sound velocity (determined from the po- 
sition of the peaks in the structure function) at the very 
same temperature. 

The very observation that A(u) — within a fi- 
nite velocity range implies a sub-exponential scaling and 
thereby suggests looking for a more accurate scaling 
Ansatz. The power-law hypothesis 

,5(2)(z,t) «t-^^(5(2)(i/r^), i<uct (4) 

proves to be the correct choice. Because of the con- 
servation law (O, we expect that temporal and spatial 
rates coincide with one another, namely 71 = 72 = 7- 
The value of 7 can be estimated from the behaviour of 
(5(2) (0, t) which appears to follow a clean power law (see 
circles in Fig. 2). A best fit (solid line) yields 0.606±0.008 
which strongly hints at 7 = 3/5. By adopting this 7 value 
in Eq. the rescaled profiles (5(2) are plotted in Fig. 3 at 
different times (without loss of generality, Q is set equal 
to 1). The scaling Ansatz (0J is finally validated by the 
extremely good overlap of the rescaled distributions in 
the interval delimited by the secondary peaks, observed 
in correspondence of the sound velocity Us- It is interest- 
ing to notice that while the system-size L = 8190 allows 
obtaining a clean scaling behavior for the perturbation 
diffusion, direct simulations of the heat conductivity in- 
dicate that L = 30, 000 is not large enough to obtain a 
looser estimate of the divergence rate 0, Q . 

In view of the conservation law expressed by Eq. Q , it 
is tempting to interpret the perturbation profile (5(2) (i,t) 
as a probability density function (pdf ) and thereby com- 
pare our results with those expected for an ensemble of 
Levy walks ^ simple schematization of a Levy 

walk consists in a particle moving ballistically in be- 
tween successive "collisions" whose time separation is 




FIG. 2: Maximum height 5(2) (0, t) of the infinitesimal pertur- 
bations profile (circles) and mean square displacement a'^{t) 
(diamonds) versus time. 

distributed according to a power law, %p{t) oc t^^^^ 
(/i > J while their velocity is chosen from some 
symmetric distribution which, in the simplest setup, re- 
duces to 'P(w) = {S{u — 1) + S{u + l))/2. The propagator 
(the pdf to find in x at time t, a particle initially localized 
at a; = 0) of such a process is [13 

F(x, t) oc < ^11 

t^-^' \x\=t 

\x\>t 

It can be easily shown that, up to the length |a;| — t, 
the propagator P{x, t) scales as in Eq. Q with the expo- 
nent 7 = The cutoff at |a;| = t is & consequence of 
the finite constant velocity and leads to sharp peaks at 
the outermost wings (see also the inset of Fig. 3, where 
the propagator for the Levy walk with jj, = I/7 — 5/3 
is compared with the direct simulation of the HPG). In 
fact, these peaks correspond to single flights during the 
observational time. The only relevant difference with the 
results of HPG simulations, namely the broad secondary 
peaks exhibited by 6(2) j disappears as soon as the 5-Dirac 
functions in the definition of 'P{u) are replaced by Gaus- 
sian distributions with a suitable width (see again the 
inset in Fig. 3). 

As for the secondary peaks, a best fit of their height 
gives a decay as i-i i5±0.03^ This is to be confronted 
with the theoretical prediction j23| which, in this 

case amounts to t~^^/^^. Much of the deviation is to be 
attributed to the not yet established asymptotic decay in 
the tail. 

Let us now discuss the connection between dynamics 
in tangent space and energy diffusion. The latter pro- 
cess can be investigated from the evolution of an initially 
localized, finite perturbation Awi(O), 

Ai?,(0 = ^(^.,(^) + A«,(0)^ - ^.,2(0 = 

m,v,(t)^v,(t) + !^(Aw,(t))2, (6) 




FIG. 3: Rescaled perturbation profiles aXt = 40, 80, 160, 320, 
640, 1280, 2560, and 3840 (the width increases with time), for 
7 = 3/5. The profiles have been obtained by averaging over 
10* realizations. In the inset, the profile At t — 640 (solid 
line) is compared with the propagators of a Levy walk for 
an exponent M = § with a fixed velocity u — 1 (dotted line) 
and a Gaussian distribution with average equal to 1 and rms 
0.036) (dashed line). In the last two cases, the propagators 
have been obtained by averaging over 10* different realiza- 
tions. 

and thereby performing an ensemble-average {AEi{t)) to 
get rid of statistical fluctuations. By directly computing 
{AEi{t)), we found that even after averaging over 10^ 
realizations, it was not possible to obtain meaningful re- 
sults on time scales shorter than those considered in the 
previous simulations. In fact, while an increasing accu- 
racy is required to resolve tinier deviations, statistical 
fluctuations remain of order one at all times. Addition- 
ally, there are fluctuation ampliflcations due to Avi sud- 
den jumpsj23]. However, in the limit of infinitesimal per- 
turbations, the probability of such jumps vanishes and 
the last term in the r.h.s. of Eq. © can be identified 
with 6(^2) ■ As a consequence, not only the total energy of 
the perturbation is connserved, but because of Eq. ||2J), 
the sums of the two terms in the r.h.s. of Eq. ^ are 
separately conserved. Finally, if the perturbation is ran- 
domly chosen, rriiViSvi) vanishes, so that the relevant 
properties of energy diffusion are captured by the pseudo- 
Euclidean norm of Svi{t). 

Therefore, from the scaling behaviour for (5(2) we can 
conclude that energy diffuses anomalously in the HPG. 
This can be also shown more directly, from the behaviour 
of the mean square displacement (T^(t) = X]j *^'^(2) (*i 
which is expected to scale as a {t) oct'^. A best fit of our 
numerical data (see diamonds in Fig. 2) yields P « 1.35, 
a value in agreement with Ref . [8b] . Since in a Levy walk, 
the exponent /3 is equal to 3— /z, we expect (3 = 4/3, which 
thus confirms the validity of the Levy walk interpretation. 
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Energy superdifFusion can be finally related to the 
anomalous behaviour of k, by exploiting the link estab- 
lished by linear response theory between heat conduc- 
tivity and the decay of spontaneous fluctuations of the 
energy density By following Ref. the diver- 

gence rate a of heat conductivity can be connected to 
the anomalous diffusivity exponent (3 through the simple 
relation a — (3 — 1 = 2 — ijl = 2 — 2/7, which implies 
a = 1/3 in our case j23 |. a value in perfect agreement 
with the prediction of dynamical rcnormalization group 
i. 

The clean results summarized in the first three figures 
of this Letter are made possible by the peculiarity of HPG 
dynamics which allows for a strong reduction of statis- 
tical fluctuations. It should not be, however, forgotten 
the absence of exponential instabilities whose presence 
would have obliged to averaging over an exponentially 
growing number of trajectories. Therefore, one might 
suspect that the Levy-walk scenario is peculiar to the 
HPG, but this hypothesis contrasts the observation that 
the HPG scaling behaviour coincides with the predictions 
based on very general hydrodynamic arguments. On the 
other hand, if Levy walks generally occur in Id systems, 
evidence for this behavior should be found in, e.g., truly 
chaotic models: however, the above mentioned numeri- 
cal difficulties strongly suggest the need to define a com- 
pletely new protocol to answer this question. 

Let us now briefly discuss the origin of the anoma- 
lous behavior. The evolution in tangent space is fully 
determined by three elements: energy plus momentum 
conservation and the correlations between collisions. In 
Ref. [23, it was introduced a Id model where the energy 
of each particle is randomly exchanged (with one of the 
two neighbors) in uncorrelated collision events. In such 
a simplified context it was rigorously proved that Fourier 
law is satisfied. It is thus natural to ask whether the 
diverging heat conductivity observed in the HPG is due 
by the additional conservation of momentum (which is 
absent in the model of Ref. 23] ). By replacing the HPG 
collision pattern with a set of completely uncorrelated 
events, we found a normal behaviour. Therefore, we are 
led to conclude that the anomaly is entirely contained in 
the long-range correlations between collisions. 

S. Lepri is acknowledged for useful discussions. This 
work is part of the project PRIN2003 Order and chaos 
in nonlinear extended systems funded by MIUR Italy. 
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